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New model equations are derived for dynamics of self- aggregation of finite-size particles. Differences 
from standard Debye-Ifiickel and Keller-Segel models are: the mobility of particles depends 
on the configuration of their neighbors and linear diffusion acts on locally-averaged particle density. 
The evolution of collapsed states in these models reduces exactly to finite-dimensional dynamics of 
interacting particle clumps. Simulations show these collapsed (clumped) states emerge from smooth 
initial conditions, even in one spatial dimension. 

Keywords: gradient flows, blow-up, chemotaxis, parabolic-elliptic system, singular solutions 



> 

o 
o 

o 
in 
o 



> 

X 
J3 



PACS numbers: 

Modeling finite size effects in the aggregation of in- 
teracting particles requires modifications of the class of 
Debye-Hiickel equations. This problem is motivated by 
recent experiments using self-assembly of nano-particles 
in the construction of nano-scale devices 0. Fundamen- 
tal principles underlying the self-assembly at nano-scales 
are non-local particle interaction and nonlinear motion 
due to variations of mobility at these scales. The model 
should account for the change of mobility due to the fi- 
nite size of particles and the nonlocal interaction among 
the particles. The local density (concentration) of par- 
ticles is denoted by p. For particles interacting pairwise 
via the potential — G(|r|), the total potential at a point r 
is $(r) = — / p(r')G(|r — r')dr' — G * p where * denotes 
convolution and G > for attracting particles. The ve- 
locity of the particle is assumed to be proportional to the 
gradient of the potential times the mobility of a particle, 
11. The mobility can be computed explicitly for a sin- 
gle particle moving in an infinite fluid. However, when 
several particles are present, especially in highly dense 
states, the mobility of a particle may be hampered by 
interactions with its neighbors. These considerations are 
confirmed, for example, by the observation that the vis- 
cosity of a dense suspension of hard spheres in water di- 
verges, when the density of spheres tends to its maximum 
value. Many authors have tried to incorporate the depen- 
dence of mobility on local density by putting n = p,{p) 
and assuming p,{p) — > as p — » pmax = 1 \4. Vanish- 
ing mobility leads to the appearance of weak solutions in 
the equations, to singularities and, in general, to massive 
complications and difficulties in both theoretical analysis 
and computational simulations of the equations. 

Alternatively, we suggest that the mobility fi should 
depend on an averaged density p over some sampling 



volume, rather than on either the potential, or the exact 
value of the density at a point. This assumption makes 
sense from the viewpoints of both physics and mathe- 
matics. From the physical point of view, the mobility of 
a finite-size particle must depend on the configuration of 
particles in its vicinity. While attempts have been made 
to approximate this dependence by using derivatives of 
the local density, this approach may lead to unphysi- 
cal negative diffusivity 5*. Hence, we assume instead 
that the local mobility depends on an integral quantity p 
which is computed from the density a&'p — H * p. Here 
H[y) is a local filter function, which may be of much 
shorter range than the potential G. Several filter func- 
tions are possible, with examples H{v) being a (5-function, 
an exponential (or inverse-Helmholtz in ID) or a top- 
hat filter. Alternatively, for crystal growth applications, 
one may assume a filter function H{v) with zero aver- 
age, so the mobility does not depend on the magnitude 
of the density, only on the relative position of particles. 
The mathematical analysis and the reduction property 
derived in this paper holds, regardless of the shape of 
the filter function H and the potential G is, so long as 
they are both nice (e.g., piecewise smooth) functions. 

Aim of the paper This paper treats the following 
continuity equation for the evolution of density, 

^ = -divJ with J = - DVp- p{-p)pV<^ . (1) 

Here Z) > is a constant, while p and <E> are defined by 
two convolutions involving, respectively, the filter func- 
tion H and the potential G, 

p ~ H * p and $ = G * p . (2) 

While model H1I2|I preserves the sign of p, we shall see 
that it allows the formation of (5-function singularities 
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even when D > and in the case of one dimension. The 
generahty and predictive power of model H1I2|I can be 
demonstrated by enumerating a few of its subcases. For 
example, when H{r) = S{r) and G{x) — e"'""'/' this sys- 
tem reduces to the generalized chemotaxis equation 
For the choice G = 5'{r), H = S'{r) one obtains a variant 
of the inviscid Villain model for MBE evolution . If the 
range of G (denoted by is sufficiently small, one may 
approximate (at least formally) the integral operator in 
<i> = G * p as a differential operator acting on the density 
p, namely, $ = p + IV ■ [IVp). As was demonstrated 
in 1^ , equation (QJ then becomes a generalization of the 
viscous Cahn-Hilliard equation, describing aggregation of 
domains of different alloys. In the case that p is constant, 
H{x ~ y) — 5{x — y) , and G is the Poisson kernel, then 
equation JQ) becomes the Poisson-Smoluchowski equa- 
tion for the interaction of gravitationally attracting par- 
ticles under Brownian motion . 

All four subcases of model (|1I2|I described above re- 
quire a singular choice of the smoothing kernels G and 
H. However, the generalized functions required in these 
singular choices may be approximated with arbitrary ac- 
curacy by using sequences of nice (for example, piecewise 
smooth) functions. We shall concentrate on cases where 
the functions H and G remain nice, and derive the results 
in this more general, more regular, setting. Thus, regu- 
larization introduces additional generality, which enables 
analytical progress and makes numerical solution of the 
equations easier. 

Gradient flow motivation Model (|1I2(I may also 
be motivated from the viewpoint of gradient flows and 
thermodynamics. The associated free energy decreases 
monotonically in time and remains finite even when the 
solutions tend to a set of delta- functions. This energy 
remains finite for all choices of p{])). In contrast, the 
traditional approach in which the mobility is assumed 
to depend on the un-smoothed density p{p) will fail by 
allowing the energy to become infinite. This additional 
energetic argument reinforces the choice of the regular- 
ized model in pi2|l . The calculation will be performed in 
arbitrary spatial dimensions. For the energy to be well- 
defined, we assume that the function G, describing the 
interaction between two particles, is everywhere positive, 
so the interaction is always attractive. In one dimension, 
we also assume that the kernel function G is symmet- 
ric and in n dimensions that the interaction is central. 
These assumptions are physically viable for all the pre- 
vious subcases. Our derivation for the free energy will 
remain valid for an arbitrary filter function H . 

We start with equation ^ for mass conservation, in 
the case when the mobility p{j}) in the particle flux J is 
not constant. We shall derive a variant of equation 
as a gradient flow defined by. 
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where {f , ^) = J d^x is the pairing, for a suitable 
test function ^E* and where E is the following energy. 



E[p]= Djp{logp^l)dx+^Jp^dx, 



(4) 



with <I> = G * p and "p = H * p. The first term in energy 
E[p] in Q decreases monotonically in time for linear dif- 
fusion of p. The second term defines the H^^ norm of 
p in ID for the choice G{x) = exp(— |a;|/a), and it pro- 
vides a generalization of this norm for an arbitrary (but 
positive and symmetric) function G. Energy E[p\ in 
remains finite, even when the solution for the density p 
concentrates into a set of (5- functions. 

The variation of the free energy E[p] in equation is 



5E[p] 



{DH *{logp) + <S>)Spd' 



where the density variation 5p is chosen in the form, 

5p = -div(p^(p)V*) , 

determined by the function 'J, which is assumed to be 
smooth, |lO|. After integrating by parts, the correspond- 
ing variational derivative is 



div(p/x(p)V((DH * (logp) + $)) ) , * 



/SE ^ 

The resulting variant of equation is, cf. (0) 

^ = - div J™ , J„, = - pp{p){DH*^Vp+V^) , (5) 
ot ^ p ' 

in which the linear diffusion of density p is mollified by 
smoothing with H . Finally obtaining the particle fiux J 
in equation |^ requires simplifying the first term of 3m 
in (O to ZJVp . Thus, only the contribution — /i(p)pV<I> 
to the particle fiux in equation is variational. 

Energetics The free energy E[p] in 10} decreases 
monotonically in time under the evolution equation 
A direct calculation yields. 



dE[p] 
dt 



I ^—\3m?d^X. 
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We shall see that the resulting mass-conserving motion 
of finite-size particles governed by © leads to a type of 
'clumping' of the density p into a support set consist- 
ing of 5-functions. One may check that this monotonic 
decrease of energy persists for more general functions, in- 
cluding these weak solutions supported on (5- functions, as 
discussed below. These weak solutions are then found to 
spontaneously emerge in numerical simulations of equa- 
tion and dominate the initial value problem. 

Dynamics of clumpons We consider motion under 
equations (|1I2I) in one spatial dimension. Substituting 
the following singular solution ansatz 
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p{x,t) = y^^w^{t)5{x ~ qi{t)) 



(7) 



i=l 
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into the one-dimensional version of either equation 
or (O with D = and integrating the result against a 
smooth test function cj) yields a closed set of equations 
for the parameters Wi{t) and qi(t), i — 1,2, N, of the 
solution ansatz Q. Namely, 
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«^z(t) = 0, ftW = -I]"'iA^(p)G'fe-<z,), (8) 



where p = Y^m=i^'mH(qj(t) - qm{t)). Thus, the den- 
sity weights Wi{t) — Wi{0) — wi are preserved for 
i — 1,2, . . . ,N , and the corresponding positions qi{t) fol- 
low the characteristics of the velocity u — ~iiij))VG * p 
along the Lagrangian trajectories defined by a; = qi{t). 
This result holds in any number of dimensions, modulo 
changes to allow singular solutions supported along mov- 
ing curves in 2D and moving surfaces in 3D. When £> 7^ 
and divDVp ~ p, the weights Wi decay essentially expo- 
nentially in time, under evolution by equation 

In Fig. 1, we demonstrate that the solutions © do in- 
deed appear spontaneously in a numerical simulation of 
equation with D = 0.02. Hence, they are essential in 
understanding its dynamics. The simulation started with 
a smooth (Gaussian) initial condition for density p. Al- 
most immediately, one observes the formation of several 
singular clumpons, which evolve to collapse eventually 
into a single clumpon. Observe that the mass of each 
individual clumpon remains almost exactly constant in 
the simulations, as required by equation with D = 0. 
Also note that the masses of two individual clumpons 
add up when they collide and "clump" together. Con- 
sequently, all the mass eventually becomes concentrated 
into a single clumpon, whose mass (amplitude) is the to- 
tal mass of the initial condition. This simulation shows 
that for a general system the long-term evolution 
evolves into the collective motion of a finite number of in- 
dividual clumpons. This conclusion is supported by anal- 
ysis showing the superposition of clumpons represented 
by the singular solution lO is an invariant manifold of 
the gradient flow equations Q|. 

Rate of blow up Analysis of the evolution of a density 
maximum reveals the clumping process results from the 
nonlinear instability of the gradient flow in equation 



, one may 



when D = 0. For a particular case G — e"'^'/" 
show that for large enough peak, a density maximum 
Pm(t) = p{xm(t),t) becomes infinite in finite time. The 
motion of the maximum is governed by 

^Pm ^ {pIi - Pm^ix,n)) > ^ " PmM) , (9) 

where M = J pdx is total mass and we have used the 
fact that satisfies $ — a^^xx = P- The last inequality 
holds, because G < 1 is bounded and p is everywhere 
positive. Thus, if at any point the maximum of p exceeds 
the (scaled) value of the total mass, then the value of 
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FIG. 1: Numerical simulation demonstrates the emergence 
of particle clumps, showing formation of density peaks in a 
simulation of the initial value problem for equation us- 
ing smooth initial conditions for density with I = 1, G{x) = 
H{x) = exp( — |a;|), p{'p) = 1. The vertical coordinate repre- 
sents J) = H * p, which remains finite even when the density 
forms 5-functions. Insert: The behavior of the inverse ampli- 
tude of the maximum peak l/pm{t) (crosses) versus analytical 
solution of equation Q (solid line). 



the density maximum Pmit) must diverge in finite time, 
which produces (5-functions in finite time. From the 
density amplitude must diverge as pm — oP' — t). To 
illustrate this divergence, we have plotted the comparison 
between the predicted collapse of 1/pm and numerics in 
the insert of Fig. 1. The formation of singularities in 
Fig. 1 occurs both at the maximum, and elsewhere, The 
subsidiary peaks eventually collapse with the main peak. 

Competing length scales of H &z G An interesting 
limiting case arises, when the scale of non-locality of H 
is much shorter than the range of the potential G. For- 
mally, this limit corresponds to H(x) — > S{x). In prac- 
tice, we may select a sequence of piecewise smooth func- 
tions H^{x) = exp(— |a;|/e)/(2e) which converge weakly to 
a (5-function. For each function H^{x) in this sequence, 
no matter how small (but positive) the value of e, the ex- 
act ODE reduction ||HJ) still holds. What happens in the 
limit of very small epsilon, as e — > 0? In investigating 
the limit as e — > 0, we performed a sequence of numer- 
ical simulations with the function G{x) = exp(— jxj/a) 
being held fixed and H^{x) = exp(— a;/e)/(2e) varying, 
for a sequence of decreasing values of the ratio e/a. This 
simulation is illustrated in Fig. 2 for e/a — 1/10, which 
demonstrates the formation of flat clumps of solutions. 
The mechanism of this phenomenon is the following. The 
vanishing mobility at p = 1 caps the maximal density at 
p = 1 in the long-term. This leads to the appearance of 
flat mesas in 'p{x) for large t. On the other hand, when 
H{x) S{x), 'p(x,t) — i- p(x,t) pointwise in x, which 
forces p{x, t) to develop a flat mesa, or plateau, structure 
in which the maximum is very close to unity, as well. This 
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is precisely what is predicted for the chemotaxis equation 
with H{x) = S{x) and n{p) = 1 - p 




FIG. 2: Evolution of a Gaussian initial condition for p(x,0) 
with /i(p) — 1 — p and H — exp( — |a::|/e)/(2e) where e — a/10. 
The solution quickly forms a plateau of maximal possible den- 
sity (pmax = !)• 

In the limit H ^ S, model (|1I2|I recovers ordinary dif- 
fusion of local density. This is a singular limit, because it 
increases the order of the differentiation in the equation. 
Since ordinary diffusion is known to prevent collapse in 
one dimension this singular limit should be of con- 
siderable interest for further analysis. 

Conclusion and Open Problems A new model was 
proposed and analyzed for the collective aggregation of 
finite-size particles driven by the force of mutual attrac- 
tion. Starting from smooth initial conditons, the solu- 
tion for the particle density in this model was found 
to collapse into a set of delta- functions (clumps), and 
the evolution equations for the dynamics of these clumps 
was computed analytically. The energy derived for this 
model is well defined even when density is supported on 
(5-functions. The mechanism for the formation of these 
(5-function clumps is the nonlinear instability governed 
by the Ricatti equation (j^, which causes the magnitude 
of any density maximum to grow without bound in finite 
time. At first sight, it may seem that the emergence of 
i5-function peaks in the solution might be undesirable and 
perhaps should be avoided. However, these ^-functions 
may be understood as clumps of matter, and the model 
guarantees that any solution eventually ends up as a set 
of these clumps. Consequently, further collective mo- 
tion of these clumps may be predicted using a (rather 
small-dimensional) system of ODEs, rather than dealing 
with the full non-local PDEs. The question of how many 
clumps arise from a given initial condition remains to be 
considered. One may conjecture that clump formation 
is extensive; so that each clump forms from the material 
within the range of the potential $, determined by G. 



On a longer time scale, the clumps themselves continue 
to aggregate, as determined by the collective dynamics 
(jHJ of weak solutions This clump dynamics is also a 
gradient flow; so that eventually only one clump remains. 

An analogy to point vortex solutions of Euler's equa- 
tions for two-dimensional ideal hydrodynamics may be 
drawn here. Prediction of motion of inviscid fluid is gov- 
erned by a set of nonlinear PDE in which pressure in- 
troduces non-locality. A drastic simplification of motion 
occurs, when all the vorticity is concentrated in delta- 
functions (point vortices) |l2l | . The motion of point vor- 
tices lies on a singular invariant manifold: if started with 
a set of point vortices, the fluid structure will remain a 
set of point vortices. However, a smooth initial condition 
for vorticity does not split into point vortices under the 
Euler motion. In the present model, though, the physical 
attraction drives any initial conditions towards a set of 
delta- functions, so one is guaranteed to obtain effectively 
finite-dimensional behavior in the system after a rather 
short initial time. 

Because two scales are present in the smoothing func- 
tions H and G, the effects of boundary conditions war- 
rant further study. For example, the H ^ S limit should 
also allow formation of boundary layers. Here, these 
boundary issues were avoided by using periodic boundary 
conditions. However, it is natural in physical situations 
to apply, e.g., Neumann boundary conditions to the par- 
ticle flux J. Thus, the issue of boundary conditions may 
become important in a more physically realistic setting. 
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